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We consider the hybrid problem of reconstructing the isotropic electric conductivity 
of a body Q from interior Current Density Imaging data obtainable using MRI mea- 
surements. We only require knowledge of the magnitude \ J\ of one current generated 
by a given voltage / on the boundary <90. As previously shown, the corresponding 
voltage potential u in f! is a minimizer of the weighted least gradient problem 



argmin{ / a(x)|Vu| : u € H 1 ^), u\gn = f}, 



> 

with a(x) = \J(x)\. In this paper we present an alternating split Bregman algorithm for 
ON . treating such least gradient problems, for a € L 2 (0) non-negative and / € -ff 1//2 (<9f2). 

We give a detailed convergence proof by focusing to a large extent on the dual problem. 
This leads naturally to the alternating split Bregman algorithm. The dual problem also 
turns out to yield a novel method to recover the full vector field J from knowledge of its 
magnitude, and of the voltage / on the boundary. We then present several numerical 
experiments that illustrate the convergence behavior of the proposed algorithm. 



h 1 Introduction 



The classical Electrical Impedance Tomography (EIT) problem seeks to obtain quantitative 
information on the electrical conductivity a of a body from multiple measurements of volt- 
ages and corresponding currents at its surface. The extensive study of this problem has 
led to major mathematical advances on uniqueness and reconstruction methods for Inverse 
Problems with boundary data. See the excellent reviews [6], jl], and [13]. However, by 
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now it is well understood that the problem is exponentially ill-posed, yielding images of low 
resolution away from the boundary [J5], [T9~] . 

A new class of Inverse Problems seeks to significantly improve both the quantitative accuracy 
and the resolution of traditional Inverse Boundary Value Problems by using data which can 
be determined in the interior of the object. These have been dubbed Hybrid Methods, as 
they usually involve the combination of two different kinds of physical measurements, and 
several recent advances are presented in the present issue of the journal. 
In this paper we continue our study of the Current Density Impedance Imaging (CDII) 
problem of reconstructing the conductivity of a body based on measurement of currents in 
its interior. Such measurements have been possible since the early 1990 due to the pioneering 
work of M. Joy's group at the University of Toronto [IB] , [T7] . The idea was to use Magnetic 
Resonance Imaging (MRI) in a novel way, to determine the magnetic flux density B induced 
by an applied current. It is important to note that in the problem addressed in this paper, 
we only require knowledge of the magnitude \ J\ of one current generated by a given voltage / 
on the boundary dQ. The analytic and numerical methods presented here do not necessarily 
depend on MRI. Since the results only require knowledge of the magnitude of one current, 
they may lead to simpler physical methodologies to obtain such data (see for instance |¥Tj). 
The problem of recovering the isotropic conductivity a of an object from knowledge of the 
magnitude of one current density \J\ in the interior was studied in [221 [231 [221 [22]- See 
[26] for a review, and for numerous references to other hybrid approaches to conductivity 
imaging. In this paper we present an alternating split Bregman algorithm for the numerical 
solution of this problem, along with a convergence proof. 

Let a be the the isotropic conductivity of an object Q C MJ 1 , n > 2 and let J be the 
current density vector field generated by imposing a given boundary voltage /. Then the 
corresponding voltage potential v satisfies the elliptic equation 

V-(«tVv) = 0, v\ 9n = f. (1) 

By Ohm's law J = —a\/v. Hence the voltage potential v satisfies the degenerate elliptic 
equation 

V- (p^ Vw ) =0. v|an = /- (2) 

In general there is no uniqueness for viscosity solutions \9\ of the above equation [361 123] • 
However, in [24J and [22] authors proved that if the potential v satisfies eqauation (1) above 
then it minimizes the energy functional E(v) = j Q |J||Vf | associated to the equation (j2J) 
and, moreover that this functional has a unique minimizer in if x (f2). Thus, the conductivity 
is uniquely determined by the magnitude of the current generated by imposing a given 
boundary voltage and the corresponding voltage potential is the unique solution of the 
(infinite-dimensional) minimization problem 

min{ / |J||Vt;| : v G H\Q). v\ dn = /}. (3) 
Jn 

A simple iterative procedure to solve (E]) was given in [23]. This method was only defined 
when Dirichlet problems such as (1) for successive approximations of a were guaranteed to 
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produce solutions with non-vanishing gradients in ft. For planar regions, this led to the 
requirement that the given boundary voltage be almost two-to-one. (See [21] for the precise 
definition) . 

In this paper we present a convergent alternating split Bregman algorithm to find the min- 
imizer of ([3]) for any given boundary voltage / G H l ^ 2 {dVt). More generally, let ft be a 
bounded region in R n , n > 2. Also let a G L 2 (ft) be a non-negative function and let 
/ G H l l 2 (dVl). Consider the minimization problem 



min / a(x)|Vv|, (4) 

{vem{n),v\an=f} Jn 

and assume that it has an optimal solution in H l (Vt). The algorithm we present in this 
paper will converge to a minimizer offjlj). 

The problem (J3J) belongs to a general class of problems of the form 

mm G(u) + F(Lu), (5) 
ueHi 

where L : Hi — > H2 is a bounded linear operator, the functions G : Hi — > R U {00} and 
F : H 2 — S'lU{oo} are proper, convex and lower semi- continuous, and Hi and H 2 are real 
Hilbert spaces. To see this, let Hi = H%(Q), H 2 = (L 2 (fi)) n , and Lu = Vw; fix u f G if 1 (ft) 
with u f \n = f and define F : (L 2 (tt)) n -»■ R and G : flj(fi) -> R as follows 



F(d) := / a|rf + Vu/|dx, G(tt) = 0. 
in 



(6) 



One approach to the problem (jSJ) is to write it as a constrained minimization problem 



min G(u) + subject to Lu = d, (7) 

ueHi,deH 2 

which leads to the unconstrained problem: 

min G (u) + F(d) + hLu-d\\ 2 . (8) 
To solve the above problem, Goldstein and Osher [12J introduced the split Bregman method: 

(u k+ \d k+1 ) = argmm u&HM {G(u) + F(d) + ±\\b k + Lu-d\\ 2 2 }, (9) 

b k+l = b k + Ly k+l _ d k+\ 

Since the joint minimization problem ([§]) in both u and d could sometimes be hard or 
expensive to solve exactly, Goldstein and Osher [T2] proposed the following alternating split 
Bregman algorithm to solve the problem (J5J) 

u k+1 = a rgmm ueHi {G(u) + ^\\b k + Lu k+1 -d\\ 2 2 }, (10) 

d k+1 = ^gmm d£H2 {F{d) + ^\\b k + Lu k+1 -d\\ 2 2 } ) (11) 

b k+l = b k + Lu k+l -d k+1 . (12) 
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Cai, Osher, and Shen [8] and independently Setzer [3U [35] proved convergence results for the 
above algorithm, under the assumption that Hi and H 2 are finite dimensional. Motivated 
by the infinite dimensional problem flU), in j2T] the authors recently presented general con- 
vergence results for the alternating split Bregman algorithm in infinite dimensional Hilbert 
spaces. In this paper we will study the following alternating split Bregman algorithm for the 
Dirichlet problem fll]). 

Algorithm 1 (Alternating split Bregman algorithm for weighted least gradient 
Dirichlet problems) 

Let u f G H\Q) with u f \ dU = f and initialize b°,d° G (L 2 (Q)) n . For k > 1: 



1. Solve 

Au k+1 = V ■ {d k {x) - b k (x)), u k+l \ dn = 0. 

2. Compute 

dk+1 = { maxllVn^ 1 + Vu } + b k \ - ^Q} |g^±|g±^ - Vu f if \Vu k +\x) + Vu/ + b k (x)\ ± 0, 
\ -Vu f if \Vu k+1 (x) + Vu f + b k (x)\ =0. 

3. Let 

b k+1 (x) = b k {x) + Vii fc+1 (x) - d k+1 {x). 



The following theorem, which we will prove in Section 3, guarantees convergence of the 
Algorithm 1. 

Theorem 1.1 Let a G L 2 (f2) be a non-negative function and f G H l l 2 (dVt). Then for 
any Uf G and any b°,d° G (L 2 (Q)) n the sequences {b k }k£N, {d k }k&N, and {u k }k£N 

produced by Algorithm 1 converge weakly to some b, d, and u. Moreover {Vu k+1 — d k } kGN 
converges strongly to zero and 

oo 

II Vu k+1 -d k f< oo. (13) 

k=0 

Furthermore u := u + Uf is a solution of the minimization problem (TJP, d = Vu, V • b = ; 
and 

b = ^(a^-) in tt\ {x : \Vu(x)\ = 0, a(x) ^ 0}. (14) 
A \ \vu\ J 

Note that convergence of the Algorithm 1 does not require uniqueness of the minimizers of 
the problem (j4]). In the case of the hybrid inverse problem where a(x) is the magnitude of 
the current density vector field J for some unknown conductivity a we can say more. 
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Theorem 1.2 Let Q be a bounded region in C M. n with connected C 1,a bounday. Assume that 
a = | J\ > a.e. in ft, where J G (L 2 (Q)) n is the current density vector field generated by an 
unknown conductivity a G C a (Q) by imposing the voltage f on dQ. Then the corresponding 
voltage potential u is the unique minimizer of (0J , and the sequences b k , d k + Vuf, andu k +Uf 
produced by Algorithm 1 converge weakly to —J/X, Vu, and u, respectively. 

Remark 1.3 A generalization of the above theorem holds when Q contains perfectly conduct- 
ing Uq or insulating Uj inclusions 12^ . The sequences b k , d k + Vw/ ; and u k + Uf produced 
by Algorithm 1 will converge weakly in Q\Uc to — J/X, Vu, and u, respectively. 

Our approach to §5§ is to start by working on the dual problem, which is well suited to a 
Douglas- Rachford splitting algorithm. This leads naturally to the alternating split Bregman 
algorithm. Indeed Setzer ([311 ES]) proved the equivalence of the two methods. In the 
process, we discovered that, for the inverse conductivity problem of interest here, the dual 
problem has a unique solution, which is in fact (up to a constant) the current J ! (See 
12. 3p . This shows that the alternating split Bregman algorithm is quite natural for the hybrid 
problem of Current Density Impedance Imaging. 

The paper is organized as follows. In Section 2 we study the dual problem. Based on the 
analysis of the dual problem, in Section 3 we present our convergence results and prove 
Theorems 11.11 and 11.21 In Section 4 we will prove that the Algorithm 1 is stable with respect 
to possible errors in solving the Poisson equations fl35|) at each step. In section 5, we present 
several numerical experiments to demonstrate the convergence behaviour of our algorithm, 
in particular when the boundary function is not two-to-one. 



2 The dual problem 

Rockafellar-Fenchel duality [31] will be the starting point for our proof of Theorem 11.11 
In this section we study the relation between the problem (j3J) and its dual problem. To 
begin fix uj G H l (Vt) with u\q = f and let a be a non- negative function in L 2 (Q). Define 
F : (L 2 {Vt)) n — )■ R and G : ->• R as follows 

F(d) := / a\d + Vu f \dx, G(u) = 0. (15) 
Jq 

Then the problem (jlj) can be written as 

(P) min G(u) + F(Vu). 

By Fenchel duality [31], the dual problem corresponding to the problem (P) can be written 



as 



(D) - min { G *(-V-b) + F*(b)}, 
be(L 2 (n)) n 

Recall that the Legendre-Fenchel transform F* of a functional F on a Hilbert space H is the 
function on H* given by 

F*(b) = sup{(d,6) -F(d) : deH}. 
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One can easily compute G* : H — > R: 

[ oo it u ^ 0. 
The following lemma provides a formula for F* : (L 2 (Q)) n — >■ R. 
Lemma 2.1 Lei a 6 L 2 (f2) fre non-negative and Uf G H 1 ^) uffl't/i = /. T/ien 

F*(b) = { ~( Vu f^ \K x )\< a ( x ) a.e. in Q, 

oo otherwise. 

Proof. Assume 

\b(x)\ > a(x), 

on a subset U of Q with positive Lebesgue measure. Then 
F*(b) = sup (d,b) - / a|d + V«/| 

= — (6, Vm/)+ sup I (d, fe) — / a\d\dx 
de(L 2 (n)) n \ Jn J 

> -(b, Vuf) + sup A / (|o| 2 - a{x)\b\)dx = oo, 
AeR Jc/ 



where for the last inequality we choose d(x) = Xb(x) for x G U, and d(x) = otherwise. 
Now assume 

\b(x)\ < a(x), a.e. 

then 



F*(b) = -(b,Vu f )+ sup \(d,b)- / a|d|d:c) 

= — (b, Vuf) + sup / (6 • d — a\d\)dx 
de(L 2 (n))™ Jn 

< — (b, Vm/) + sup / \d(x)\(\b(x)\ — a(x))dx 

de(L 2 (n))"- Jn 

< ~(b,Vu f ). 
Finally note that taking d = in (fT7|) yields 

F*(6) > -(6,V«/). 

It follows from ( I16p that the dual problem can be explicitly written as 

max{< Vu f ,b>: b G (L 2 (fi)) n , |6(ar)| < o(ar) a.e. and V ■ b = 0}. 
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(17) 



□ 



Now let u G iZj(ft), b G {L 2 {Vt)) n with V • b = 0, and \b\ < a a.e. in tt. Then 

/ a\Vu + Vu f \dx > / \b\\Vu + Vu f \dx > / b ■ (Vu + Vu f )dx =< b,Vu f > . (19) 
Jn Jn Jn 

Hence if we denote the optimal values of the primal and dual problem by v (P) and v(D), 
respectively, then v(P) > v(D). This is a general fact. Moreover, since both of the functions 
F and G are continuous, it follows from Theorem 4.1 in [TT] that strong duality holds, i.e., 
v(P) = v(D) and the dual problem (D) has an optimal solution. 

The algorithm we propose seeks to construct a solution of the dual problem (D). In the rest 
of this section we show how this leads to the solution of the primal problem (P) as well. 



Vit/ + {p ■ p € (L 2 (f2)) n , p = mb on {x : a(x) > 0},m G M^} if \b\ < a a.e. 

otherwise, 



Lemma 2.2 Let F* be defined as IfTb)) . Then 
dF*(b) = 
where 

Mf, := {m : Q, — > [0, oo) [ m is measurable and m(x) = if \b(x)\ < a(x)}. 

Proof. Assume \b\ < a a.e. in Q. Let m G and consider d G (L 2 (Q)) n with |d| < a a.e 
in Q. We have 

(mb,d) < / m\d\\b\dx = / m\d\\b\dx = / m\b\ 2 dx = / m\b\ 2 dx = (mb, 6), 

Jn J{\b\=a} J{\b\=a} Jn 

and hence 

F*(d) - F*(b) = -(Vu f , d) + (Vu /( 6) > (-Vu/ +mb,d-b). 

Note that the points where a(x) = (hence also b(x) = d(x) = 0) do not contribute to any 
of the terms above. On the other hand 

F*(d) - F*(b) > (-Vu f + mb,d-b), 

trivially holds if \d\ > a on a set of positive measure, as F*(d) = oo in that case. Therefore 

-Vu f + mb G dF*(b). 

Now let p G dF*(b) and define p := Vu^ + p. Then for any d G (L 2 (f2)) n with |d| < a a.e in 
F*(d) - F*(b) = -(Vu f , d) + (Vu f , b) >{p,d-b). 

Consequently 

(p,d)<(p,b), (20) 
for all d G (L 2 (Q)) n with < a(x). In particular if we let 

l!i| ^ 

3 otherwise, 
7 



then it follows from ( 120|) that 

(p,b) > (p,d) = I a\p\dx > I |6||p|dar> / b-p>(p,b). (21) 



Therefore all inequalities in ( 12 ip are equalities. Hence there exists a non-negative function 
m such that 

if \b(x)\ < a{x) 
mb if \b{x)\ = a(x) ^ 0. 



p{x) = 

This completes the proof. □ 



Now we are ready to prove the following proposition which is a special case of the Rockafellar- 
Fenchel duality theorem [31] in convex analysis. Given one solution of the dual problem (D), 
this result gives a description of all the solutions of the primal problem (P). 

Proposition 2.1 Let b be an optimal solution of the dual problem (D). Then u G Hq(Q) is 
an optimal solution of the primal problem (P) if and only if 

VuedF*(b). (22) 

Proof. Let u be a solution of the primal problem. Then, as we saw in ffTOl) 

/ a\Vu + Vu f \dx > / \b\\Vu + Vu f \dx 
Jn Jn 

> b- (Vu + Vu f )dx = (b,Vu f ). 
Jn 

Since the duality gap is zero, both inequalities are equalities and 

b(x) = a(g) -||±|gL if | Vu + Vtt,| ^ 0. (23) 

Therefore for x with a(x) ^ 0, 

Vu{x) = —Vuf + m(x)b(x), 

where 

m(x) = < a ( x ) 1 ;| r 

I otherwise. 

In view of Lemma [2.21 we conclude Vu G dF*(b). 

Now assume V« G F*(b). Then by Lemma 12.21 there exists m G M~ h such that for x with 
a(x) 7^ 0, Vu(x) = — Vuf + m(x)b(x). Therefore 

(b,Vu f ) = / b ■ (V« + Vu f )dx = / \b\\Vu + Vu f \dx = / a\ Vu + Vu f \dx, 
Jn Jn Jn 
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which means u is a minimizer of the primal problem (P). 



□ 



In the next section, we shall see that the above relation ( 122"]) between optimal solutions of 
the primal and dual problems is at the heart of Algorithm 1. The relation (123]) shows that 
b is determined on the set where \Vu(x) + Vu/(x)\ does not vanish. We record this fact as 
a separate partial uniqueness result: 

Proposition 2.2 Let u be an optimal solution of the primal problem and assume b\ and b\ 
are two optimal solutions for the dual problem (D). Then 

6i = S 2 in Q\{x: \Vu{x) + Vu f {x)\ = 0,a(x) ^ 0}. (24) 

Proof. By Lemma I2TT1 Vtt G dF*(bi) C\dF*(bi). hence it follows from lemma I2T21 that there 
exist m% G and m 2 G Mg such that 

Vu = — Vuf + mi&i and Vii = — Vuj + 
Thus mibi = m^- Since b\ and Si are both optimal, by ( 12 3 j) 

| Si (a;) | = |S a (ar)| = a(x) on \ {x : |Vu(x) + V«/(ar)| = 0, a(x) 7^ 0}. 

Hence (El]) follows. □ 
If a(x) is the magnitude of the current corresponding to a conductivity a the above propo- 
sition yields uniqueness of solutions to the dual problem. 

Corollary 2.3 Let Q be a bounded region in C IR n with connected C 1,a bounday. Assume 
a = I J I > a.e. in fi, where J G (C Q (f2)) n is the current density vector field generated by an 
unknown conductivity a G C a (fl) by imposing the voltage f on dQ. Then the corresponding 
voltage potential u is the unique minimizer of (j3J) and the dual problem (D) has a unique 
solution b. Furthermore b = — J. 

Proof. The proof follows from Proposition 12.21 . the uniqueness Theorem 1.3 in ([23]) and 
equation (T23l . □ 

Remark 2.4 The above corollary generalizes to the case when Q contains perfectly con- 
ducting (Up) and/or perfectly insulating Uj inclusions. Proposition \2.S\ together with the 
uniqueness result in 12^ show that a solution b of the dual problem (D) equals the current J 
in Q\U p . We omit the details. 



3 Convergence analysis 

In this section we present proofs of Theorems 11.11 and 11.21 The proofs rely on the repre- 
sentation of solutions of the primal problem (122]) in Proposition 12.11 Notice that the dual 
problem (D) can be written in the form of an inclusion problem 

0eA(b) + B(b), (25) 
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where A := <9(G*o(— V*)) and B := OF* are maximal monotone operators on (L 2 (Q)) n . If 
we can compute a solution b of the problem (1251) as well as d G B(b) = dF*(b) then, by 
Proposition 12. 1} u = V _1 (<i) will be a solution of the primal problem (P). The Douglas- 
Rachford splitting method in Convex Analysis yields precisely such a pair (b,d). Following 
this route to the primal problem leads naturally to the Alternating Split Bregman algorithm, 
as we explain below. 

Let if be a real Hilbert space and let A, B : H — > 2 H be two maximal monotone (set valued) 
operators. For a set valued function P : H — > 2 H , let Jp denote its resolvent i.e., 

J P = (Id + P)- 1 . 

The sub-gradient of a convex, proper, lower semi-continuous function is maximal monotone 
and if P is maximal monotone then Jp is single valued [TJ [31] . Lions and Mercier [T8] 
showed that for any general maximal monotone operators A, B and any initial element Xq 
the sequence defined by the Douglas-Rachford recursion: 

x k+1 = (J a (2Jb ~ Id) +Id- J B )x kl (26) 

converges weakly to some point x G H such that p = Jp{x) solves the inclusion problem 
( 125]) . Recently Svaiter [38] proved that the sequence p k = J V B(xk) also converges weakly 
to p. This fact will be important for our problem. The following theorem describes the 
Douglas-Rachford splitting algorithm and summarizes these known convergence results. 

Theorem 3.1 Let H be a Hilbert space and let A, B : H — ^ 2 H be maximal monotone 
operators and assume that a solution of (25)) exists. Then, for any initial elements xo and 
Po and any X > 0, the sequences p k and x k generated by the following algorithm 

= JxA&Pk -x k ) +x k - Pk 
Pk+i = J\B(xk+i), (27) 

converge weakly to some x and p respectively. Furthermore, p = J\p,(x) and p satisfies 

G A(p) + B(p). 

We wish to apply the Douglas-Rachford splitting algorithm to the operators A := d(G*o(— V*)) 
and B := dF*. We need an efficient way to evaluate the resolvents J\A{2pk ~ x k) an d 
J\B{xk+i) at each iteration. If we let 

x k = X(b k + d k ), p k = \b\ fc>0, (28) 

then to evaluate J\A(^Pk~ x k) an d J\B{%k+i) f° r & h k > 0, it is enough to find the minimizers 
u h+1 and d k+1 of the functionals 

i?(u) =|| Vu + b k -d k || 2 , (29) 

and 

I*(d) = [ a\d + Vu f \dx+ || b k + Vu k+1 - d || 2 (30) 
Jn 
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and set b k+1 = b k + Vu k+1 — d k+1 . Indeed the resolvents J\A^Pk ~ x k) an d J\B(%k+i) can 
be computed as follows 

J XA (2 Pk - x k ) = \{b k + Vu k+l - d k ), 

and 

JxB(x k+1 ) = \(b k + Vu k+1 -d k+1 ), 

(see [34], [35] for a proof). Finding the minimizer of (l29|) amounts to solving a Poisson equation 

Au k+1 =\7-(d k - b k ), u k+1 \ an = 0. 
As well, the minimizer of the functional I k (d) can be computed explicitly as follows 

dk+ i. = { max{| V^ +1 + Vu f + b k \ - f , 0} ~ Vu f if \Vu k +\x) + Vu f + b k (x)\ 0, 

-Vu f if \Vu k+1 (x) + Vu f + b k (x)\ =0. 

We are thus led to the Algorithm 1 for simultaneously finding solutions of both the primal 
problem and the dual problem. To prove the strong convergence of the series (JTBl we will 
need the following simple lemma on firmly non-expansive operators. 

It is well known that the operator T := Ja{2Jb — Id) + Id — Jb is firmly non-expansive, i.e., 
T = \ld + ~R such that R is non-expansive: 

|| Rx — Ry ||<|| x — y \\ for all x,y G H. 

Lemma 3.2 If T : H — )■ H is a firmly non-expansive operator and x k +\ = T(x k ) with 
Xq G H , then 

II %k+i ~ x || 2 + || Xk+l - x k \\ 2 <\\ x k - x || 2 . 

Proof. Since T is firmly non-expansive, R = IT — Id is a non-expansive operator. 
Hence 

|| Rxk — Rx || 2 = — || x k — x || 2 +2 || Txk — Tx || 2 —2 || (Id — T)xk — (Id — T)x || 2 . 
Therefore we have 

i(|| Xk - x || 2 - || Rxk -Rxr)=\\x k -xr- || x k+1 - x || 2 - || x k+l - Xk || 2 . 

Since R is non-expansive, the left hand side of the above inequality is non-negative, which 
completes the proof. □ 

Proof of Theorem 11.11 By interpreting Algorithm 1 as a Douglas- Rachford splitting algo- 
rithm as detailed above, weak convergence of the sequences d h , and b k follows immediately 
from Theorem 13. 11 To prove the estimate (fT3"j) . let T = J\aC^J\b — Id) + Id — J\b- Since T 
is firmly non-expansive, by Lemma [3.21 we have 

|| x k+ i - x || 2 + || x k+ i - x n || 2 <|| x k - x || 2 , (31) 
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where x is the weak limit of x k with T(x) = x. By the above inequality we have 



^ II %k+i - x k || 2 < oo. (32) 

fc=0 

Now observe that 

x k+1 - x k = X(b k+l + d k+l - b k - d k ) = X(Vu k+1 - d k ), 
and hence ( ITBl) follows. 

By Theorem 13.11 p = Xb is a minimizer of the dual problem and J\QF*{X(d + b)) = Xb. 
Therefore 

Xb + XdF*(Xb) = X(d + b) & dedF*(Xb). 

In view of (fT3l) the sequence {u k } k eN is bounded in Hq(£1) and therefore it has a weakly 
converging subsequence. Let u be a weak cluster point of the sequence {u k }. Then Vw = 
d G df*(p) and hence by Proposition 12.11 w is a solution of the primal problem (P). On the 
other hand Vu = d for every weak cluster point u of the sequence {u k } k ^N- Since V is 
injective on Hq{VL), {u k } ke ^ has at most one weak cluster point. In view of Proposition I2.2[ 
the proof is now complete. □ 

Proof of Theorem 11.21 The proof follows from ll.ll combined with the uniqueness Theorem 
1.3 in [21], and our Corollary E33 □ 



4 Approximate alternating split Bregman algorithm 

In this section we show that the alternating split Bregman algorithm converges to the correct 
solutions even in the presence of possible errors at each step in solving Poisson equations. 
The proof relies on the following theorem about the Douglas-Rachford splitting algorithm. 

Theorem 4.1 (Svaiter |38j) Let X > 0, and let {a k } k£ N and {(3 k } k eN be sequences in a 
Hilbert space H . Suppose G ran(A + B), and ^fcejvdl ak II + II Pk II) < 00 ■ Take xq G H 
and set 

x k+ i = x k + J lA {2{J XB x k + n ) - x k ) + a k - (J XB Xk + 0k), k > 1. (33) 

Then x k and p k = J\BX k converge weakly to x G H and p G H , respectively and p = J\bx G 
(A + B)-\0). 

The proof of the above theorem in infinite-dimensional Hilbert spaces is due to Svaiter [5B] 
(see also [7]). It suggests the following approximate version of our algorithm 

Approximate alternating split Bregman algorithm: 

Let u f G H l (Q) with u f \ 9U = f and initialize b°,d° G (L 2 (fi))". For k > 1 
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1. Find an approximate solution u k of 

Au k+1 = V ■ {d k {x) - b k {x)), u k+1 \ dn = 0, 
with || Vu k — Vu k x || < ctfc, where u^ 1 is the exact solution of the above problem. 

2. Compute 



d 

3. Set 



, = f maxHV^ 1 + Vu f + b k \ - f , 0} ^tv^l ~ Vu / if |V^ +1 (*) + Vu, + 6 fc (x)| ^ 0, 
\ -Vtt/ if \X7u k+1 (x) + Vu f + b k (x)\ = 0. 



By Theorem 14. II and an argument similar to that of Theorem II. II we can prove the following 
theorem about convergence of the sequences u k , d k , and b k produced by the above algorithm 

Theorem 4.2 Let a e L 2 (VL) be a non-negative function and f E H 1 ^ 2 ^). If 

oo 

^2 a k < OO, 

k=l 

then for any Uf G H 1 ^) and any b°,d° G (L 2 (Q)) n the sequences {b k }k£N, {d k }kGN, and 
{u k } k^N produced by the approximate alternating split Bregman algorithm converge weakly 
to some b, d, and u. Moreover {Vu k+1 — d k }k^N converges strongly to zero and 

oo 

|| Vu k+l -d k || 2 < oo. 

k=0 

Furthermore u := u + Uf is a solution of the minimization problem $4$, b is a solution of the 
dual problem (D), d = Vm ; V • b = 0, and 

S = Ifa-^-] in Q\{x: |V«(ar)| = 0, a(x) ^ 0}. (34) 



There is also an analogue to Theorem 11.21 for the Approximate alternating split Bregman 
algorithm, which we omit. 



5 Numerical experiments 

In this section we study numerically the convergence behavior of the proposed alternating 
split Bregman algorithm. In a model problem, we cseek to reconstruct the conductivity 
from knowledge of the magnitude of one current density |J| given inside the unit square 
Q = (0, 1) x (0, 1) and the corresponding voltage potential / on dQ. In [21] authors presented 
a simple iterative algorithm to recover the conductivity a from the knowledge of (| J\, f). To 
be well defined, this algorithm requires the solution u and all intermediate functions u k 
produced by the algorithm to have non-vanishing gradients in Q. The alternating split 
Bregman algorithm proposed in this paper does not require |Vm| > in Q and converges 
(by Theorem 11.11) to an optimal solution of (jl]). We perform several numerical experiments 
to illustrate this convergence behavior. 
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5.1 Data simulation 

To simulate the internal data \J\, we use an abdominal human CT image rescaled to a 
realistic range of tissue conductivity, with values varying from 1 to 1.8 S/m The scaled 
conductivity distribution, on a uniform grid 128 x 128, is shown in Figure [TJ 



Figure 1: The original conductivity distribution used in the data simulation. 

Given the above conductivity distribution, we first solve numerically the Dirichlet problem 

V-aVv = 0, v\dtt = f, 

on the grid 128 x 128. To provide high accuracy of the numerical solution, we look for a 
solution of the form v = Uh + u, where Uh is the harmonic function satisfying the Dirichlet 
boundary condition and u is the solution to the Poisson equation 



with the zero boundary conditions. For the range of conductivity values described above, 
the norm ||w|| is small compared to that of \\uh\\- Such a representation is also helpful in 
the Algorithm 1; the function Uh needs to be computed just once. As forward solvers, we 
use the FORTRAN sofware FISHPACK and MUDPACK for elliptic problems (see [37]). We 
run those routines with the double precision. Comparison of the numerical solutions with 
an analytical one (in the case of constant conductivity) verified that the relative L 2 -error 
does not exceed 10~ 6 . Note that if such an error were to exceed 10~ 5 , then it may severely 
affect the quality of the reconstructed images. We combined the above solvers with the 
numerical differentiation via the three- or five-point Lagrangian interpolation, to preserve 
the high accuracy needed in the reconstruction algorithms. Once the solution u is computed 
the magnitude of the current density in Q, which is the data for our problem, is simulated 
as \ J\ = er|Vw|. 




— V ■ aVu = V ■ crWuh 
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5.2 Reconstruction algorithms 



For reconstruction we use the alternating split Bregman algorithm, Algorithm 1, choosing for 
Uf the harmonic extension Uh of the boundary voltage /. This leads to a simpler algorithm. 
For comparison purpose, we will also apply the so-called simple iterations algorithm from 
[23]. Below we outline both of algorithms for convenience. 

Algorithm 2 (Simplified Alternating split Bregman algorithm) 

Let Uh be the harmonic extension of / to Q and initialize d°, b° G (L 2 (Q)) n . 

1. Beginning with k — 0, solve 



Au = V ■ (d k (x) — b k (x)), u\ dn = 



(35) 



and let v 



U + Uh- 



2. Compute 




Vv k+1 (x)+b k (x) 
|Vu fc + 1 (x)+6*(a;)| 



if \\7v k+l (x) +v k (x)\ ^ 0, 
if \Vv k+l (x) + b k (x)\ = 0. 



3. Let 




b k (x) +Vv k+1 {x) -d k+l {x). 



J 



a = 



|Vv fc+ ir 



Otherwise, := A; + 1 and repeat the process. 



The simple iterations algorithm. 



Let Uh be the harmonic extension of / in Vt and initialize uq = uj,a± 



\J\ 



1. Beginning with k = 1, solve the problems 



Am = V • (a k S/u h ), u\ d n = 0. 



V k = U + Uh- 



2. Update 



J 



Cfe+l — 



|Vu fc r 
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3. If -gjfcU < Tol, then STOP. Otherwise, let k := k + 1 and repeat the process. 

In both algorithms, the harmonic part Uh of the solution is computed via the over-relaxation 
method, and the Dirichlet problem for the Poisson equation is solved numerically by using 
the implicit conjugate- gradient method (see [23]) in which the iterative process is constructed 
by minimizing the error of each approximation in the energy norm. The correction vector 
from the Krylov space is determined on each iteration. By virtue of the implicit method, a 
five-diagonal matrix is inverted on each iteration using a preconditioner. 



5.3 Numerical reconstructions 

We use A = 1 in all the numerical experiments with the alternating Bregman algorithm. 
In our first experiment we choose the almost two-to-one boundary voltage f(x,y) = y. 
The results obtained by applying the split Bregman algorithm with N = 1, 5, 10, 30, 50, 100 
iterations are shown in Figure [2J A larger image of the conductivity reconstructed using the 
alternating split Bregman algorithm with N = 60 is shown in Figure [3J This image may be 
compared with the original image in Figure [TJ 




Figure 2: Conductivity reconstruction using the alternating split Bregman algorithm with 
N = 1, 5, 10, 30, 50, 100 iterations (shown from the left upper corner to the right lower corner) 
for the almost two-to-one boundary condition f(x,y) = y. 
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Figure 3: Conductivity reconstructed using the alternating split Bregman algorithm with 
N = 60 iterations for the almost two-to-one boundary condition f(x,y) = y. 




Figure 4: Reconstruction using the simple iterations algorithm with N = 1, 5, 10, 30, 50, 100 
iterations (show from the left upper corner to the right lower corner) for the almost two-to- 
one boundary condition / = y. 

For comparison we repeat the above experiment with the same almost two-to-one boundary 
condition using the simple iterations algorithm. Figure H] shows the resulting conductivity 
for different number of iterations. For two-to-one boundary data, the results of the two 
algorithms are similar, although the simple iterations method slightly outperforms the alter- 
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nating Bregman algorithm in this case. 



The Tables [T] and [2] give the numerical errors of experiments with alternating split Breg- 
man and simple iterations algorithms for different levels of tolerance and the computation 
times on a Dell Precision T5400 workstation with an Intel Xeon 64-bit 2 core processor. 



Tolerance 


Relative L 2 — error vs EXACT 


Elapsed time(s) 


Total number of iterates 


5xl0- b 


0.0156 


703 


122 


lxlO" 4 


0.0148 


574 


99 


2xl0~ 4 


0.0075 


443 


76 


5xl0" 4 


0.0166 


276 


47 



Table 1: Numerical errors and elapsed times for alternating split Bregman 



Tolerance 


Relative L 2 — error vs EXACT 


Elapsed time(s) 


Total number of iterates 


5xl0- 5 


0.0030 


672 


110 


LrlO" 4 


0.0030 


583 


99 


2xl0" 4 


0.0137 


446 


73 


5xl0~ 4 


0.0141 


264 


43 



Table 2: Numerical errors and elapsed times for simple iterations 



Recall that the simple iterations algorithm requires |Vw fc | > on Q for all k > 1. In general 
if the boundary data / is not two-to-one then there exist x G VL such that V«(i) = 0. This 
may lead to divergence of the simple iterations algorithm. For the next experiments, we 
chose the boundary voltage f(x,y) = y + 2sin(77n/), which is not two-to-one. As shown in 
Figure [5] the resulting surface z — \J\ touches the xj/-plane. For this boundary data, the 
simple iterations algorithm breaks down. However as shown in Figure [6] the alternating split 
Bregman algorithm converges. 
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Figure 5: Magnitude of the current density \J\ for the non two-to-one boundary data 
f(?,v) = y + 2sm(7ny). 

In Figure [7J we show the rate of convergence of the alternating split Bregman (Solid) and 
simple iterations (Dashed) algorithms for two-to-one boundary data f(x,y) = y (left) and 
non two-to-one data f(x,y) — y + 2sin(7iiy). There is no dashed curve on the left, as 
the errors rapidly exceed the scale in the figure. Our numerical experiments indicate that 
simple iterations slightly outperform the alternating split Bregman algorithm for two-to-one 
boundary data. However, for non two-to-one boundary data simple iterations may diverge, 
while the alternating split Bregman always converges (see Figure [7]). 

In [22] the more genral problem of recovering an isotropic conductivity outside some perfectly 
conducting or insulating inclusions was considered. The data was, as before, the magnitude of 
one current density field |J| in the interior of Q.. We proved that (except in some exceptional 
cases) the conductivity outside the inclusions, and the shape and position of the perfectly 
conducting and insulating inclusions are uniquely determined by the magnitude of the current 
generated by imposing a given boundary voltage. Since the relevant minimization problem 
is still of the form [3] in this case, the split Bregman algorithm can be applied. Figure [8] shows 
the conductivity constructed using 100 iterations of the alternating split Bregman algorithm 
in the presence of perfectly conducting (right) and insolating (left) inclusions. 
In additional experiments we examined the effect of noise in our reconstruction. The noise 
model we used is a simple stochastic model \J\ n = \ J\ + 7 * R, where R is the normally 
distributed pseudo-random matrix of the order as | J\ with mean zero and standard deviation 
of one, and 7 > is the model standard deviation chosen as 7 = 5 * |||J|||/||i?||, where S 
is the noise level, i.e., (| J\ n — \ J\)/\ J\- In Figure [9j we show reconstructed images obtained 
by 20 iterations of the alternating split Bregman algorithm for three different levels of noise. 
The numerical values of the Z 2 -relative errors vs the exact solution for Figure M are shown in 
TalbelEJ 
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Figure 6: Conductivities constructed using the alternating split Bregman algorithm with 
N — 1, 5, 10, 30, 50, 100 iterates (shown from the left upper corner to the right lower corner) 
for the non two-to-one boundary data f(x, y) = y + 2 sin(77n/). 



50 100 150 200 

Number of Iterates 



50 100 150 200 

Number of Iterates 



Figure 7: Rate of convergence for alternating split Bregman (Solid) and simple iterations 
(Dashed) for two-to-one boundary data f(x,y) = y (left) and non two-to-one data f(x,y) = 
y + 2sin(77ry). 



Low Noise (Level=0.01) 


Moderate Noise ( Level=0.035) 


Higher Noise ( Level=0.06) 


0.026 


0.080 


0.152 



Table 3: Numerical error for alternating split Bregman 
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Figure 8: Reconstruction in the presence of the perfectly conducting (right) and insulating 
(left) inclusions 





Figure 9: Low noise (left), moderate noise (middle), and higher noise (right). 
Conclusion 

We presented a convergent alternating split Bregman algorithm for least gradient problems 
with Dirichlet boundary data. This, in particular, leads to a convergent algorithm for re- 
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covering an isotropic conductivity a form the knowledge of the magnitude \J\ of one current 
generated by a given voltage / on the boundary. Duality plays an essential role in our con- 
vergence proof and leads to a novel method to recover the full vector field J from knowledge 
of its magnitude, and of the voltage / on the boundary. The alternating split Bregman 
algorithm presented here converges for non two-to-one boundary data as well as two-to- 
one boundary data. Several successful numerical experiments demonstrated the convergence 
behavior of the proposed algorithm. 
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